Fast 3D gravity and magnetic modelling using midpoint quadrature and 2D FFT

To avoid the problem of the traditional methods consuming large computational resources to calculate the kernel matrix and 2D discrete convolution, we present a novel approach for 3D gravity and magnetic modelling. This method combines the midpoint quadrature method with a 2D fast Fourier transform (FFT) to calculate the gravity and magnetic anomalies with arbitrary density or magnetic susceptibility distribution. In this scheme, we apply the midpoint quadrature method to calculate the volume element of the integral. Then, the convolution of the weight coefficient matrix with density or magnetization is efficiently computed via the 2D FFT. Finally, the accuracy and efficiency of the proposed algorithm are validated by using an artificial model and a real topography model. The numerical results demonstrate that the proposed algorithm’s computation time and the memory requirement are decreased by approximately two orders of magnitude compared with the space-wavenumber domain method.

Fast 3D gravity and magnetic modelling using midpoint quadrature and 2D FFT Xulong Wang 1 , Jianxin Liu 1,2,3 , Jian Li 1,2* & Hang Chen 4 To avoid the problem of the traditional methods consuming large computational resources to calculate the kernel matrix and 2D discrete convolution, we present a novel approach for 3D gravity and magnetic modelling. This method combines the midpoint quadrature method with a 2D fast Fourier transform (FFT) to calculate the gravity and magnetic anomalies with arbitrary density or magnetic susceptibility distribution. In this scheme, we apply the midpoint quadrature method to calculate the volume element of the integral. Then, the convolution of the weight coefficient matrix with density or magnetization is efficiently computed via the 2D FFT. Finally, the accuracy and efficiency of the proposed algorithm are validated by using an artificial model and a real topography model. The numerical results demonstrate that the proposed algorithm's computation time and the memory requirement are decreased by approximately two orders of magnitude compared with the space-wavenumber domain method.
The efficient and accurate forward modelling of the gravity and magnetic potential fields, generated by arbitrary mass density or magnetization distribution of the 3D geological body is the foundation of the geophysical applications 1,2 . Gravity and magnetic forward modelling, which can be solved in the space or Fourier domains, is part of the potential field data processing.
The earliest numerical modelling of the potential field is mainly studied in the space domain. For the space domain forward techniques, convolution integrals satisfied by the potential fields can be directly calculated by analytical solution 3-8 , Gaussian quadrature method 9-11 , Cauchy-type integral method [12][13][14] or multilevel fast multipole method 1,15 . However, whether performed by analytical formulea or numerical integration, several calculations are required for every subdivided volume cell 16 . Moreover, the discrete convolution of the weight coefficients with density or magnetization takes much time. Therefore, 3D gravity and magnetic numerical modelling suffer from a prohibitive computational cost, especially when the number of observation sites and source elements are huge 14,15 .
Another method to evaluate the convolution integrals is transformed into the Fourier domain. Many researchers have considered fast computation of integral kernel, including the calculation of the integral in the Fourier domain [17][18][19][20] and by discretization of the operator and calculation in the space domain [21][22][23][24][25] . Then, Li et al. 17 and Dai et al. 18 developed the space-wavenumber domain method based on 2D Gauss-FFT 26 to calculate gravity and magnetic anomalies. Their research mainly focus on reducing the boundary effect of the standard FFT in the Fourier domain, which significantly improving the accuracy.
Forward modelling for geophysics kernels has utilized the FFT in a variety of situations. Vogel 27 first applied the FFT algorithm to the integral kernel where observations and dipoles are aligned on a regularly-spaced grid. For a multilayer model, Zhang and Wong 21 introduced the use of the Block-Toeplitz-Toeplitz-Block (BTTB) structure for fast 3D gravity forward modelling and inversion. Wu 22 combined the Gaussian quadrature and FFT algorithm to develop a fast method for complex 3D gravity forward modelling. Subsequently, Chen and Liu 23 used an analytical expression to calculate the weight coefficients, but only for gravity case. Yuan et al. 25 achieved the magnetic forward modelling based BTTB matrix. Hogue et al. 24 discussed the properties of the coefficient matrix of the gravity and magnetic potential field. However, the majority of previous work apply the analytical soulution or Gaussian quadrature method to caluculate the kernel martrix which hit a computational cost and memory requirement constraint, rendering them prohibitive in scenarios with a large number of mesh elements.
In this study, we propose a novel method to solve large-scale gravity and magnetic potential field forward problems swiftly. In this scheme, we introduce an approach to explore the symmetric structures of the coefficient matrix for the gravity and magnetic kernel. To reduce the computation time, the midpoint quadrature 28,29 is used to calculate the coefficient matrix, and the 2D FFT is applied to achieve the fast discrete convolution of the coefficient matrix with the density or magnetization function. The accuracy and efficiency of our algorithm are examined using a synthetic model and a realistic topography model.

Basic theory
Statement of problem. For general mass density or magnetization distribution, the integral response of both gravity and magnetic signals can be expressed as the convolution of the Green's function G(r, r ′ ) with the source function f (r ′ ) 30 which can be further simplified as where r is the field coordinates, r ′ is the source coordinates, h(r) represents the response of the gravity or magnetic signals, J(r, r ′ ) = cG(r, r ′ ) is the weight coefficients matrix. For the gravity forward problem, Numerical integration of 3D integral. The computational burden of evaluating element integral by an analytical solution 7,31 or numerical integration 14,15 is pretty large, even if the density or magnetization profile of the cell is as simple as constant. To solve it, we introduce the midpoint quadrature method 28,29 to calculate the volume element of the integration by multiplying the value of the midpoint with the length of the integration interval. In this section, J(r, r ′ ) can be computed by applying the midpoint quadrature method. The model And, we assume that the mass density or magnetization of each grid point is constant. Therefore, for the fields point (x l , y m , z 0 ) on the 2D horizontal grid, the general 3D convolution problem can be written as is the density or magnetization of the source point. Equation (3) can be further expressed as Taking gravity field's Green function as an example, the weight coefficient of the formula (6) can be expressed as Similarly, we can obtain Green's function the weight coefficients for gravity and magnetic anoamly quantities. The final source matrix f 2L×2M is obtained by filling the original source matrix f L×M with zero element, which is equal to the size of kernel matrix.
Appling the 2D inverse Fourier transform on Eq. (9), the gravity anomaly in the space domain can be obtained where F −1 2D denotes the 2D inverse discrete Fourier transform operators. ⌊⌋ L×M means to take the first L row and M column elements of the matrix. Finally, the gravitational anomaly is obtained by accumulating the influence of the N layer source on the observation height. Repeating the above process, we can also easily get the magnetic anomaly.

Results
In this section, we examine the proposed method by presenting two different kinds of numerical experiments. At first, we use a benchmark scenario model 32,33 to verify the proposed algorithm and evaluate the performance of the proposed algorithm by comparing it with existing implementation algorithms. Then, we apply the proposed approach to a local real topographic model in New Zealand to demonstrate the practical application. All the numerical experiments are implemented on a laptop with Intel Core i5-10210U 1.6 GHz and 8 GB of RAM. Fig. 1, the benchmark model is used to verify the correctness and evaluate the performance of the gravity and magnetic forward modelling code. The model setup is identical to that used in Pan et al. 32 , with a modelling domain of 100 km × 100 km × 100 km , which contains two cubic bodies. The model region is discretized into 128 × 128 × 128 cells with an equal interval of 781.25 m in three directions. The number of observation points is 128 × 128 on a horizontal observation plane with a height of 12.5 km.

Performance comparison with existing implementations. As shown in
The forward results in comparison with the analytical solution 30 are shown in Figs. 2 and 3. For the gravity case, the maximum errors for the gravity field g z and gravity gradient T g zz are 1.07 × 10 −5 mGal and 1.06 × 10 −6 E, respectively. For the magnetic case, the maximum errors for the magnetic field B z and magnetic gradient T m www.nature.com/scientificreports/ are 1.7 × 10 −5 nT and 5.91 × 10 −6 nT/m, respectively. The trivial difference between the numerical results and analytical solution for all data types obviously verifies the accuracy of the proposed algorithm.
To evaluate the numerical performance of the proposed method, the forward results are compared with the space-wavenumber domain method 18 , which transforms the 3D integral equation methods into a 1D integral with different wavenumbers. In addition to the method mentioned 18 , we also compare our method with the compact difference scheme for solving the Posson's equation 32 .
The gravity field g z and gravity gradient T g zz of the analytical solutions are shown in the first columns of Fig. 2. The absolute differences of our method and the space-wavenumber domain method with 8 gauss nodes 18 are shown in the second and third columns of Fig. 2, respectively. For the gravity field g z and the gravity gradient T g zz , we can find that the differences of our method are two orders of magnitude lower than those of Dai et al. 18 . Taking T g zz as an example, the execution time, rrms, maximum absolute error and computing resources for the different methods are shown in Table 1. The maximum error and the relative root mean square error (Rrms) of our method is 1.06 × 10 −6 E and 7.89 × 10 −6 % , which are much lower than Dai et al. 18 (1.34 × 10 −5 E and   18 and Pan et al. 32 , respectively. Our method has a speedup ratio of 78 and 124 over the Dai et al. 18 and Pan et al. 32 methods, respectively. It can be seen that our proposed method has higher computational accuracy and efficiency. We also study the execution time and memory requirement in cases of different model sizes. The model is discretized into different small cells, which are from 50 × 50 × 50 cells to 250 × 250 × 250 cells. As illustrated in Fig. 4b, the computation time of our method increases from 0.029 s (50 × 50 × 50) to 1.81 s (250 × 250 × 250). By contrast, the method of Dai et al. 18 spends 0.9 s and 699.26 s for the same condition. It can be evident that the computation time of both approaches increases exponentially with the increase of the mesh as shown in Fig. 4b.   www.nature.com/scientificreports/ However, the computation time of our method is more gently compared with Dai et al. 18 . Moreover, our method requires only 1/16 of the computation memory of Dai et al. 18 in Fig. 4a. So, it can be further considered that our method is more favor for 3D large-scale gravity and magnetic forward problems.
Topography effect on airborne data. Figure 5 shows the final testing model, which derived from a digital elevation model (DEM) of Rotorua city, New Zealand, with significant terrain. Our motivation is to set up a rugged topography model to investigate the effect of topography on airborne gravity and magnetic data. The DEM data with a spatial resolution 3 ′′ × 3 ′′ , can be obtained from the http:// dwtkns. com/ srtm/. The DEM data are resampled with a spacing of 50 m in the horizontal directions. The area's maximum elevation of the undulating terrain is 1091 m, and the lowest is 3.2 m. Thus, the rugged topography model can be built by 400 × 400 × 100 cells. For all elements below the topographic ground surface, we assign a constant value of 2760 kg m −3 . The magnetic field strength is 53861 nT, and the magnetic declination and inclination are 20.9 • and − 63.5 • , respectively. The local magnetic data can be obtained from the Magnetic Declination website (Magnetic declination, n.a.). The magnetic susceptibility of the region is taken as 0.02 SI. We set up a horizontal computational plane with 400 × 400 observation sites that cover the entire area at a local height of 2000 m to imitate the airborne gravity and magnetic survey.   www.nature.com/scientificreports/ We calculate the vertical gravity to evaluate the accuracy and efficiency of the proposed method based on the terrain mode as shown in Fig. 5b. Figure 6 compares the vertical gravity computed using our method with the analytical solution based on cell merging and parallel computing 33 , and both solutions are in excellent agreement. The maximal relative error of the two methods is less than 2.2 × 10 −3 % . Additionally, compared to the 46,779.56 s required by the analytical solution, our method only took 2.56 s. Our solution has a speedup ratio of 18,273 over the analytical solution based on cell merging and parallel computing. The excellent agreement of the two solutions, as well as the considerable speedup clearly demonstrate the better performance of our method for tackling complicated topography model.
In addition, we also compute the gravity fields (Fig. 7a,b), the gravity gradient tensors (Fig. 7c-i), the magnetic fields ( Fig. 8a-c) and the magnetic gradient tensors (Fig. 8d-i) caused by the terrain model. From Fig. 7f, we can see that the maximum value of T g xx + T g yy + T g zz is less than 1.0 × 10 −7 E, which is approximately equal to 0. Obviously, it again demonstrates that our method is suitable for a complex topography model.

Discussion
In this article, we propose an efficient approach by utilizing the midpoint quadrature for computing the weight coefficient matrix and the 2D FFT method to directly calculate gravity and magnetic anomalies. We evaluate its precision and effectiveness using an artificial model and a real topography model. The results show that our method can almost ignore the effect near the boundary compared with the traditional Fourier domain method based on Gauss-FFT. At the same time, this approach not only greatly improves efficiency, but also requires much less computational memory. The results show the great potential of our method to be applicable to 3D fast inversion of gravity and magnetic data in a large-scale realistic terrain model. www.nature.com/scientificreports/

Data availability
The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.